%matplotlib inline
from scipy import linalg
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML, displaySciPy - scientific library in python
Adapted from the notebooks by
Introduction
SciPy builds upon NumPy.
Among others, SciPy deals with:
- Integration (scipy.integrate)
- Optimization (scipy.optimize)
- Interpolation (scipy.interpolate)
- Fourier Transform (scipy.fftpack)
- Signal Processing (scipy.signal)
- Linear Algebra (scipy.linalg)
- Sparse matrices (scipy.sparse)
- Statistics (scipy.stats)
- Image processing (scipy.ndimage)
- IO (input/output) (scipy.io)
Animation display with python
Animation with matplotlib
You can use FuncAnimation to animate a sequence of images:
%%capture
fig, ax = plt.subplots()
xdata, ydata = [], []
(ln,) = plt.plot([], [], "ro")
def init():
ax.set_xlim(0, 2 * np.pi)
ax.set_ylim(-1, 1)
return (ln,)
def update(frame):
xdata.append(frame)
ydata.append(np.sin(frame))
ln.set_data(xdata, ydata)
return (ln,)
ani = FuncAnimation(
fig,
update,
interval=50,
frames=np.linspace(0, 2 * np.pi, 100),
init_func=init,
blit=True,
)display(HTML(ani.to_html5_video()))Another version exists (using JavaScript):
%%capture
fig, ax = plt.subplots()
xdata, ydata = [], []
(ln,) = plt.plot([], [], "bo")
ani = FuncAnimation(
fig,
update,
interval=50,
frames=np.linspace(0, 2 * np.pi, 100),
init_func=init,
blit=True,
)display(HTML(ani.to_jshtml()))More information: http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/
Animation with plotly
matplotlib works fine for advanced tuning, but is harder for simple tasks. So just try plotly for basic animations:
import plotly.express as px
from plotly.offline import plot
df = px.data.gapminder()
fig = px.scatter(
df,
x="gdpPercap",
y="lifeExp",
animation_frame="year",
animation_group="country",
size="pop",
color="continent",
hover_name="country",
log_x=True,
size_max=55,
range_x=[100, 100000],
range_y=[25, 90],
)
fig.show("notebook")See more here: http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/ ## Linear algebra
scipy for linear algebra : use linalg. It includes functions for solving linear systems, eigenvalues decomposition, SVD, Gaussian elimination (LU, Cholesky), etc.
Documentation: http://docs.scipy.org/doc/scipy/reference/linalg.html
Solving linear systems:
Find x such that: A x = b for specified matrix A and vector b.
A = np.array([[1, 0, 3], [4, 5, 12], [7, 8, 9]], dtype=float)
# A = np.random.randn(500, 500)
# b = np.ones((500, 1))
b = np.array([[1, 2, 3]], dtype=np.float64).T
# b = np.array([[1, 2, 3]], dtype=np.float64).T
print(A)
print(b)
x = linalg.solve(A, b)
print(x)
print(x.shape)
print(b.shape)[[ 1. 0. 3.]
[ 4. 5. 12.]
[ 7. 8. 9.]]
[[1.]
[2.]
[3.]]
[[ 0.8 ]
[-0.4 ]
[ 0.06666667]]
(3, 1)
(3, 1)
Check the result at a given precision (different from ==)
np.allclose(A @ x, b, atol=1e-14, rtol=1e-15)True
Remark: NEVER (or you should really know why) invert a matrix. ALWAYS solve linear systems instead!
Eigen values / vectors
A v_n = \lambda_n v_n with v_n the n-th eigen vector and \lambda_n the n-th eigen value. The associated python functions are eigvals and eig:
A = np.random.randn(3, 3)
A = A + A.T
evals, evecs = linalg.eig(A)
print(evals, "\n ------\n", evecs)
# V = [v_1, ... v_n] , columns = eigen vectors
# A = V diag(s_1,...,s_n) V^T
np.allclose(A, evecs @ np.diag(evals) @ evecs.T)[-2.24628369+0.j 3.66928214+0.j 1.02170698+0.j]
------
[[ 0.33074892 0.74511509 -0.57914476]
[ 0.85959506 -0.49114261 -0.14097965]
[ 0.38948874 0.45120111 0.80294214]]
True
EXERCISE : Eigen values/vectors
Verify numerically that the output from linalg.eig are indeed approximately eigen values and eigen vectors of matrix A above.
Hint: use https://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html
If A is symmetric you should use eigvalsh (H for Hermitian) instead: this is more robust and leverages the structures (you know they are real!)
Matrix operations
linalg.trace(A)# tracelinalg.det(A)# determinantlinalg.inv(A)# Inverse, consider NEVER using it though :)
Norms
print(linalg.norm(A, ord="fro")) # fro for Frobenius
print((np.sum(A ** 2)) ** 0.5)
print(linalg.norm(A, ord=2))
print((linalg.eigvalsh(A.T @ A) ** 0.5))
print(linalg.norm(A, ord=np.inf))4.421912143758986
4.421912143758986
3.6692821424278055
[1.02170698 2.24628369 3.66928214]
4.501267540013564
EXERCISE : Norms computation
Check numerically what is the instruction linalg.norm(A, ord=np.inf) computing. Double check with the help and a numerical test.
A = np.random.randn(3, 3)
print(linalg.norm(A, ord=np.inf))3.1174332797374276
Random generation, distributions, etc.
(see more on this at https://albertcthomas.github.io/good-practices-random-number-generators/ or https://numpy.org/doc/1.18/reference/random/legacy.html#numpy.random.RandomState)
seed = 12345
rng = np.random.default_rng(seed) # can be called without a seed
rng.random()0.22733602246716966
Visualization of various popular distributions are available here (see widgets): https://github.com/josephsalmon/Random-Widgets
Optimization
Goal: find functions minima or maxima Documentation: http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html
from scipy import optimizeFinding (local!) minima
def f(x):
return 4 * x ** 3 + (x - 2) ** 2 + x ** 4
def mf(x):
return -(4 * x ** 3 + (x - 2) ** 2 + x ** 4)
xs = np.linspace(-5, 3, 100)
plt.figure()
plt.plot(xs, f(xs))
plt.show()
Default solver for minimization/maximization: fmin_bfgs (see https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm)
x_min = optimize.fmin_bfgs(f, x0=-4)
x_max = optimize.fmin_bfgs(mf, x0=-2)
x_min2 = optimize.fmin_bfgs(f, x0=2)
plt.figure()
plt.plot(xs, f(xs))
plt.plot(x_min, f(x_min), "o", markersize=10, color="orange")
plt.plot(x_min2, f(x_min2), "o", markersize=10, color="red")
plt.plot(x_max, f(x_max), "|", markersize=20)
plt.show()Optimization terminated successfully.
Current function value: -3.506641
Iterations: 7
Function evaluations: 16
Gradient evaluations: 8
Optimization terminated successfully.
Current function value: -6.201654
Iterations: 5
Function evaluations: 12
Gradient evaluations: 6
Optimization terminated successfully.
Current function value: 2.804988
Iterations: 7
Function evaluations: 16
Gradient evaluations: 8

EXERCISE : Bassin of attraction
Draw the points on the curves with two different colors : - orange for the points leading to find the left local minima - red for the points leading to the right local minima.
Find the zeros of a function
Find x such that f(x) = 0, with fsolve.
omega_c = 3.0
def f(omega):
return np.tan(2 * np.pi * omega) - omega_c / omega
x = np.linspace(1e-8, 3.2, 1000)
y = f(x)
# remove vertical lines when the function flips sign
mask = np.where(np.abs(y) > 50)
x[mask] = y[mask] = np.nan
plt.plot(x, y)
plt.plot([0, 3.3], [0, 0], "k")
plt.ylim(-5, 5)
optimize.fsolve(f, 0.72)
optimize.fsolve(f, 1.1)
optimize.fsolve(f, np.linspace(0.001, 3, 20))
np.unique(np.round(optimize.fsolve(f, np.linspace(0.2, 3, 20)), 3))
my_zeros = (
np.unique((optimize.fsolve(f, np.linspace(0.2, 3, 20)) * 1000).astype(int)) / 1000.0
)
plt.figure()
plt.plot(x, y, label="$f$")
plt.plot([0, 3.3], [0, 0], "k")
plt.plot(my_zeros, np.zeros(my_zeros.shape), "o", label="$x : f(x)=0$")
plt.legend()
plt.show()

Parameters estimation
from scipy.optimize import curve_fit
def f(x, a, b, c):
"""f(x) = a exp(-bx) + c."""
return a * np.exp(-b * x) + c
x = np.linspace(0, 4, 50)
y = f(x, 2.5, 1.3, 0.5) # true signal
yn = y + 0.2 * np.random.randn(len(x)) # noisy added
plt.figure()
plt.plot(x, yn, ".")
plt.plot(x, y, "k", label="$f$")
plt.legend()
plt.show()
(a, b, c), _ = curve_fit(f, x, yn)
print(a, b, c)
2.5053828094575428 1.3837737554917904 0.5296873565282768
Displaying
plt.figure()
plt.plot(x, yn, ".", label="data")
plt.plot(x, y, "k", label="True $f$")
plt.plot(x, f(x, a, b, c), "--k", label="Estimated $\hat{f}$")
plt.legend()
plt.show()
Rem: for polynomial fitting, one can directly use numpy
x = np.linspace(0, 1, 10)
y = np.sin(x * np.pi / 2.0)
line = np.polyfit(x, y, deg=10)
plt.figure()
plt.plot(x, y, ".", label="data")
plt.plot(x, np.polyval(line, x), "k--", label="polynomial approximation")
plt.legend()
plt.show()/home/jsalmon/anaconda3/envs/peerannot/lib/python3.10/site-packages/IPython/core/interactiveshell.py:3505: RankWarning:
Polyfit may be poorly conditioned

Interpolation
from scipy.interpolate import interp1d, CubicSpline
def f(x):
return np.sin(x)
n = np.arange(0, 10)
x = np.linspace(0, 9, 100)
y_meas = f(n) + 0.1 * np.random.randn(len(n)) # add noise
y_real = f(x)
linear_interpolation = interp1d(n, y_meas)
y_interp1 = linear_interpolation(x)
cubic_interpolation = CubicSpline(n, y_meas)
y_interp2 = cubic_interpolation(x)
plt.figure()
plt.plot(n, y_meas, "bs", label="noisy data")
plt.plot(x, y_real, "k", lw=2, label="true function")
plt.plot(x, y_interp1, "r", label="linear interp")
plt.plot(x, y_interp2, "g", label="CubicSpline interp")
plt.legend(loc=3)
plt.show()
Images
from scipy import ndimage, misc
img = misc.face()
type(img), img.dtype, img.ndim, img.shape
print(2 ** 8) # uint8-> code sur 256 niveau.
n_1, n_2, n_3 = img.shape
np.unique(img)
# True image
plt.figure()
plt.imshow(img)
plt.axis("off")
plt.show()/tmp/ipykernel_56976/4081896913.py:3: DeprecationWarning:
scipy.misc.face has been deprecated in SciPy v1.10.0; and will be completely removed in SciPy v1.12.0. Dataset methods have moved into the scipy.datasets module. Use scipy.datasets.face instead.
256

RGB decomposition
fig, ax = plt.subplots(3, 2)
fig.set_size_inches(9, 6.5)
n_1, n_2, n_3 = img.shape
ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256))
ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256))
ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256))
plt.tight_layout()
print(img.flags) # cannot edit...
img_compressed = img.copy()
img_compressed.setflags(write=1)
print(img_compressed.flags) # can edit now
img_compressed[img_compressed < 70] = 50
img_compressed[(img_compressed >= 70) & (img_compressed < 110)] = 100
img_compressed[(img_compressed >= 110) & (img_compressed < 180)] = 150
img_compressed[(img_compressed >= 180)] = 200
plt.figure()
plt.imshow(img_compressed, cmap=plt.cm.gray)
plt.axis("off")
plt.show() C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : False
WRITEABLE : False
ALIGNED : True
WRITEBACKIFCOPY : False
C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : True
WRITEABLE : True
ALIGNED : True
WRITEBACKIFCOPY : False

Convert a color image in grayscale
plt.figure()
plt.imshow(np.mean(img, axis=2), cmap=plt.cm.gray)
plt.show()
Blurring (fr: floutage)
img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.imshow(img_flou, cmap=plt.cm.gray)
ax.axis("off")
plt.show()Widget on blurring bandwidth (cannot be rendered online for the moment, only locally)
%matplotlib widget
import ipywidgets as widgets
# set up plot
img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.axis("off")
plt.show()
@widgets.interact(sigma=(0.1, 200, 0.1))
def update(sigma=2):
"""Remove old lines from plot and plot new one"""
# [l.remove() for l in ax.lines]
img_flou = ndimage.gaussian_filter(img, sigma)
ax.imshow(img_flou, cmap=plt.cm.gray)Changing colors in an image
import pooch
import os
url = "https://upload.wikimedia.org/wikipedia/en/thumb/0/05/Flag_of_Brazil.svg/486px-Flag_of_Brazil.svg.png"
name_img =pooch.retrieve(url, known_hash=None)
img = (255 * plt.imread(name_img)).astype(int)
img = img.copy()
plt.figure()
plt.imshow(img[:, :, 2], cmap=plt.cm.gray)
fig, ax = plt.subplots(3, 2)
fig.set_size_inches(9, 6.5)
n_1, n_2, n_3 = img.shape
ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256), density=True)
ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256), density=True)
ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256), density=True)
plt.tight_layout()

EXERCISE : Make the Brazilian italianer
Create a version of the Brazilian flag as follows:

More on colors/ images:
http://josephsalmon.eu/enseignement/Montpellier/HLMA310/matplotlib_slides.pdf
Additional lectures
- http://www.scipy.org - The official web page for the SciPy project.
- http://docs.scipy.org/doc/scipy/reference/tutorial/index.html - A tutorial on how to get started using SciPy.
- http://scipy-lectures.github.io - Scipy lectures
- https://github.com/scipy/scipy/ - The SciPy source code.